Quantum thermal transport from classical molecular dynamics 
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Using a generalized Langevin equation of motion, quantum ballistic thermal transport is obtained 
from classical molecular dynamics. This is possible because the heat baths are represented by ran- 
dom noises obeying quantum Bose-Einstein statistics. The numerical method gives asymptotically 
exact results in both the low-temperature ballistic transport regime and high-temperature strongly 
nonlinear classical regime. The method can be thought of as a semi-classical approximation to the 
quantum transport problem. A one-dimensional quartic on-site model is used to demonstrate the 
crossover from ballistic to diffusive thermal transport. 

PACS numbers: 05.60.Gg, 44.10.+i, 63.22.+m, 65.80.-|-n 



Many approaches have been used to study lattice heat 
transport in bulk materials and nanostructures. For bulk 
materials, the standard method is that of Peierls based 
on Boltzmann equation for phonons [l|, 0]. For quasi- 
one-dimensional systems and nano junctions, a variety of 
techniques has been used, such as molecular dynamics 
(MD) 1^, 01, mode-coupling theory H,nonequilibrium 
Green's function (NEGF) method 0, i, H 0, [H, 
Schrodinger e quat ion method quantum Langevin 

dynamics [l^, Tl3 | , rigorous Boltzmann equations , 
etc. One of the outstanding problems in heat transport is 
to reconcile the ballistic nature at low temperatures and 
diffusive transport at high temperatures. As far as we 
know, the methods mentioned above work only in either 
ballistic regime, or diffusive regime, but none correctly 
in both. 

Molecular dynamics has the potential to be such a uni- 
versal method for heat transport. However, since MD is 
based on classical Newtonian mechanics, the quantum ef- 
fect is completely absent. Thus we can not expect that 
it is still correct at low temperatures. In fact, due to 
very high Debye temperatures for carbon based mate- 
rials, even 300 K is considered a low temperature. The 
kinetic theory of heat transport for phonons gives a for- 
mula for the thermal conductivity as k = ^cvl, where c 
is heat capacity, v is sound velocity, and I is mean free 
path. The reduction of thermal conductivity at low tem- 
peratures is mainly due to much reduced quantum heat 
capacity c, but a classical MD can only produce a con- 
stant heat capacity. 

Can we simulate a quantum system within MD? At 
first sight, this seems impossible, since classical dynam- 
ics can only produce classical results. In this paper, we 
show that the heat transport problem in junction sys- 
tems can be studied with a classical generalized Langevin 
dynamics using a quantum heat bath derived from Bose- 
Einstein statistics. Instead of the generic Nose-Hoover 
heat bath, it is essential to use the generalized Langevin 
dynamics with memory kernel and colored noises to take 
care correctly the effect of the baths. The heat baths are 



modeled as infinite numbers of coupled harmonic oscil- 
lators. A remarkable feature of the proposed dynamics 
is that it reproduces the quantum ballistic results at low 
temperature when nonlinearity can be neglected, as well 
as gives a correct high-temperature, strongly nonlinear 
result. This appears to be the only method that is nu- 
merically exact in both limits. Although the classical 
and quantum generalized Langevin equations are well- 
known, it is somewhat surprising that they are seldom 
used in molecular dynamics. In fact, they have much 
better properties with respect to heat baths; we advocate 
their use for thermal transport problems. Our method 
is inspired by the NEGF approach [13] to heat transport 
and also the quantum Langevin approach 1J| to the same 
problem. 

In the rest of the paper, we first introduce the model 
and give the equations involved. We then compare the 
MD results with Landauer formula and with the non- 
linear NEGF results. We treat a one-dimensional (ID) 
quartic nonlinear onsite model, in which we have seen 
ballistic transport at temperatures below 200 K, and dif- 
fusive transport about 1000 K for lattice sizes up to 4096. 

The general setup of our system consists of a central 
junction region connected to two semi-infinite harmonic 
lattices which serve as leads. The Hamiltonian of the 
system is 

H=Y.H^ + {u'^fv'^''u^ + {u^fV^V + H,,, (1) 

a=L,CM 

where Ha = ^{u")'^u" + ^{u")'^ K^u" , u" is a column 
vector consisting of all the displacement variables in re- 
gion a {— L, C, R), and is the corresponding conju- 
gate momentum. The superscript T stands for matrix 
transpose. We have chosen a renormalized displacement 
where nij is the mass associated with j-th 
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degree of freedom, Xj is the actual displacement having 
the dimension of length. K" is the spring constant ma- 
trix and V^'~' ~ (V^*^^)"^ is the coupling matrix of the 
left lead to the central region; similarly for V'^^. The 
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equations of motions are of the form 



il^ = 


-K^u^ + 




(2) 


ii^ = 






(3) 


zi^ = 


-K^'u^ - 




(4) 



The heat-bath degrees of freedom and can be ehm- 
inated by solving them in terms of the central variables 
and initial conditions, given, e.g., for the left lead: 



u\t) 



g{t,t')V^^u^{t')dt' 



to 



^J^u\t,)-g{tM)u\t,), (5) 
Ota 



where g{t, t') is the time-domain retarded surface Green's 
function of the left lead obtained by the solution of 



+ g{t,t')K^ ^ -5{t~t')I, 



(6) 



with the condition g{t, t') = if t - i' < 0. 

Substituting the formal solutions of the leads into 
the central region, we obtain the following generalized 
Langevin equation 14 , 3, l3] for the central part of the 
degrees of freedom: 

= -K^u^ + F„{u^-f t')u^(t')dt' + a + Cfl, 

J to 

(7) 

where F„ is the nonlinear force, S is the retarded self- 
energy of the leads, S = S/, -|- Sjj, as used in the NEGF 
calculation, but in the time domain; Sz, = V'~^^gV^'~^ . 
A similar equation holds for the right lead S/j using the 
right lead surface Green's functions. Contribution from 
the left lead due to the initial conditions is 



The expression for the right lead is analogous. The 
initial time to will be set to — oo. Using the concept of 
adiabatic switch-on, at time — oo, the three subsystems, 
left lead, central region, and right lead, are decoupled and 
the leads are in respective thermal equilibrium. We turn 
Eq. ([7]) into a stochastic differential equation by requiring 
that u^{to) and u^{tQ) are random variables. 

So far we have treated the system as a classical sys- 
tem. However, at this point, we'll make a departure and 
treat the leads quantum-mechanically. Since the lead 
system is linear, the classical equation of motion and 
quantum Heisenberg equation of motion are identical. 
At time to — oo, the leads are isolated. We assume 
that the leads obey a quantum Bose-Einstein statistics. 
This induces a random variable i^l(^) having zero mean, 
= 0; the following correlation matrix 

i^LmUff) = V^'^(^g{t,to){u^{to)u^{tof)g{t',tof 



-g{t,to){u''{to)u^{tof)g{t',tof 
-git,to){u''{to)u'^itof)g{t\tof 

+ g{t,to){u'^ito)uitof)g{t\tof)v''''.{9) 

For a sensible heat bath, the correlation should be time 
translationally invariant and independent of tQ. Indeed, 
great simplification can be done if we use the eigenmode 
representation for the matrix g: 

g{tX)^s^g's, gf = ^eit-t' f''''^^'-''\ (10) 

where S is the orthogonal matrix that diagonalizes K^, 
SK^S^ = ri^, il^ is a diagonal matrix with diagonal ele- 
ments LOjS are the positive eigen frequencies. Substi- 
tuting this result into the correlation expression, also us- 
ing the quantum equilibrium correlation values for {uu^) , 
(iiu^), (uii^), and (uu^), we obtain 



{amL{t'y) = v'^'^s'Dsv 



rLC 



(11) 



where D is a diagonal matrix with elements Dj = 
(2/K) + l)25-cosa;j(t-t')+2T^sin^,(t-t')- /M - 
[exp(/3/,?iw) — l] is the Bose-Einstein distribution func- 
tion at the temperature of the left lead. 

The appearance of the imaginary number in the last 
term in D seems ominous, as we cannot simulate a heat 
bath with imaginary correlation. The imaginary part 
comes from the fact that in quantum mechanics, ^{t) 
and ^(t') are non-commuting, and the product of the 
two is not a Hermitian operator. Such a difhculty can 
be easily overcome if we use a symmetrized correlation 
+ This amounts to interchang- 

ing t and t' and taking the transpose. The final effect is 
simply to drop the imaginary term. 

Then the question arises that such a treatment will 
not give correctly the quantum results. It turns out that 
it causes no problem, at least for the expression of heat 
current. We can show rigorously that, with the sym- 
metrized heat baths, we reproduce exactly the Landauer 
result with Caroli- formula as the transmission coefficient. 
However, the symmetrization does have a consequence to 
the quantum heat-current fluctuations. 

Using the (surface) density of states, the expression 
can be further simplified to get a rather compact result 
for the spectrum of the noises 14| . 
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FM = J {^Lml{0))e'^'dt = {^hiu;} + ^JnrM 

(12) 

where TlIco] = ^(SlM - ElM^) = -2lmV^^g[u]V^^ . 
The spectrum function F[uj] is even in lu and is a sym- 
metric matrix. Classical limit is obtained if we take 
(/(w) -t- l/2)?i w /csTi/w, where fc^ is the Boltzmann 
constant and is the temperature of left lead. 
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The thermal current in steady state can be computed 
in several equivalent ways: 



It 



-In 



dt 

= ~{u^{tfB{t))^{u^\trB{t)), (13) 



u 



where B{t) = - Y.L(t,t')u^ {t')dt' + ClW- 

The stochastic differential equation, Eq. ([7]), can be 
solved numerically in a straightforward way. Both the 
memory function (retarded self-energy S) and noise spec- 
trum F can be obtained through the surface Green's 
function g. Efficient recursive algorithms exist for the 
solution of 5 13, A set of past coordinates, u'-^(t), 



needs to be stored, in order to perform a numerical in- 
tegration due to the self-energy. We can use a simple 
rectangular rule for the integration. The random noises 
can be generated using a spectrum method Let 
the discrete Fourier transform of ^(t) be rjk = 
ttk + ibk, k = -M/2,- • - ,-1,0, 1, . . .,M/2- 1; where M 
is the number of sampling points in the discrete Fourier 
transform. Then the noises can be generated by tak- 
ing real numbers and bk {k > 0) as independent 
Gaussian random numbers with zero mean and vari- 
ance ■^F[ujk]h'M , where h is the integration step size, 
LUk = 2nk/{hM). The noise values at required times 
are obtained by an inverse fast Fourier transform as 
^(i = hi) = j^^i^rikGxp{—i2nlk/M). The numerical 
integration of Eq. ([7]) is not substantially more expensive 
than standard MD. This is because the forces are usually 
short-ranged; we only need to do the extra work for these 
sites that are directly connected to the leads. Note that 
the matrix elements of V^'" and V^"^^ are zero except 
those that have a direct connection between the center 
and leads. The computational complexity becomes more 
favorable as the system becomes larger. 

To illustrate the general method, we consider a simple 
ID model with a quartic on-site potential {(f>'^ model). 
Such a model is known to have diffusive transport in the 
classical limit [20|. The equation of motion is given by 



Kuj^i - {2K + Ko)uj + Ku 



(14) 



where the nonlinear term is nonzero only in the central 
region, i.e., /ij = /i if 1 < j < iV and /ij — otherwise. 
The required surface Green's function can be obtained 
analytically in frequency domain as g[u}\ = —A/if, where 
A is the root of the quadratic equation, K\~^ -f- (w + 
i0+)2 -2K - Kq + KX = Q, such that |A| < 1. We use 
the following expression for the heat current 



rMD 



Y\m +Uj+i)(Uj 



))• 



(15) 



Due to energy conservation along the chain, one can show 
that Eq. (|15p is equal to that defined by Eq. in steady 
state. 




FIG. 1: Thermal conductance a for the ID onsite model 
without the nonUnear interaction, with spring constant K — 
1.0eV/(amuA^), A'o = O.IK. The smooth curve is the Lan- 
_ dauer formula result, while the symbols are MD results with 



a size N = 8. The time-step h = 
steps each are used. 



10-"s and 5 x 10** MD 



We now present our numerical results. First, when 
there is no nonlinear interaction, fij = 0, the heat current 
can be computed exactly through the Landauer/Caroli 
formula, II = ^ J^^^ dLohLoTviG''TLG-Ti^){fL - /r), 

where C = {G")^ = {{to + iO+f - - E)"\ The 
molecular dynamics with the quantum heat bath repro- 
duces this result exactly. In Fig. [Tl we present the com- 
parison of MD and the exact curve. The conductance is 
defined by 



lim 



Tr 



(16) 



A numerical finite-difference with 10-percent above or be- 
low the average temperature is used. Within MD sta- 
tistical errors (computed from statistical fluctuations of 
multiple runs), the agreement is perfect. For the ballistic 
transport, the thermal conductance is independent of the 
lengths N of the system. 

A nontrivial result is obtained when the system has 
nonlinear interactions. This is presented in Fig. [21 With 
a nonlinear strength of /i = 1 [eV/(amu^A^)], we ob- 
tain quantitatively correct picture of ballistic transport 
at low temperatures and small sizes (a oc iV") and dif- 
fusive transport at high temperatures and large sizes 
(cr oc 1/iV). The low-temperature results can be com- 
pared with the NEGF ones. This is presented as smooth 
curves in Fig. [2l The NEGF results are obtained with a 
mean- field approximation to the self-energies [lo| . The 
Green's functions are iterated in equilibrium and the con- 
ductance is calculated with an approximate formula for 
the transmission, f[uj] = iTr{G''(rL + lTn)G''TR} + 
^Tt{G''TlG'^{Tr + ir,,)}, where the nonlinear effect 
is reflected in the extra nonlinear self-energy, r„ = 
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T(K) 

FIG. 2: Thermal conductance a for the ID onsite model with 
the nonlinear interaction /i = 1 eV/ (amu'^A*), spring constant 
K = 1.0eV/(amuA2), Ko = G.IK. The smooth curves are the 
NEGF results for sizes A'^ = 4 and 32, respectively, while the 
symbols are MD results with size N from 4, 16, 64, 256, 1024, 
to 4096, from top to bottom. The time-step h = 10~^^s and 
10* MD steps each are used. 



because there is very little difference between a quan- 
tum and classical system if the system is linear. The 
dynamics is such that it smoothly crosses over to the 
classical regime. Thus the method produces correctly 
results both in the quantum ballistic limit and classical 
diffusive limit. We have applied the method to a simple 
ID onsite model. Clearly, it is of general applicability. 
For example, we can use the approach to study ballis- 
tic and diffusive thermal transport in carbon nanotubes 
and graphene ribbons. We can also study the nonlinear 
effect in interfaces. The present method opens new way 
for studying quantum transport and nonlinearity. 

The author thanks Jingtao Lii, Jian Wang and Imam 
Makhfudz for discussions. The computations were per- 
formed on the clusters of the Center for Computational 
Science and Engineering and of Singapore-MIT Alliance, 
as well as on IBM cluster of Institute of High Perfor- 
mance Computing. This work is supported in part by 
a Faculty Research Grant of the National University of 
Singapore. 



~ S^). The MD and NEGF resuhs agree with each 
other at the low-temperature side very well. Clearly, the 
nonlinear NEGF results are not exact at high temper- 
atures. Thus the deviation between MD and NEGF is 
understandable. If classical heat baths are used, then 
as the temperature decreases, the thermal conductance 
increases monotonically to a size-independent ballistic 
value of (w,„ax-Wmin)fcB/(27r), where w^ax - t^min is the 
phonon band width. At the intermediate range of tem- 
peratures, no reliable methods exist that can be com- 
pared with the quantum MD results. Thus, in this dif- 
ficult temperature range, the MD results are the only 
numbers to offer. Whether we see ballistic or diffusive 
transport in a given temperature is determined by the 
mean free path of the phonons in comparison with the 
system size N . From the data in Fig. [21 we can judge 
that the mean free path is about 10'^ lattice spacings in 
temperature range of 1000 K. 

The dynamics also gives correctly the quantum average 
energy and quantum heat capacity (say, with equal tem- 
peratures for the two leads). This is consistent with the 
fact that quantum conductance is calculated correctly. In 
classical simulation, the average kinetic energy gives the 
local temperature of the system, (u|) = fc^T. However, 
this is not true in our dynamics and the kinetic energy 
is several times larger than implied by the equipartition 
theorem. Interestingly, in the limit of high temperatures 
of several thousand Kelvin, the equipartition theorem is 
restored. 

In summary, we showed that a generalized Langevin 
dynamics as a classical stochastic differential equation 
can reproduce quantum ballistic transport if the heat 
baths follow the quantum prescription. This is achievable 
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